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We present preliminar results and tests of a new general relativistic code to simulate the 
hydrodynamic collapse of a 21 solar masses star. We have assumed spherical symmetry 
and used the formalism of Misner and Sharp to construct a finite-difference scheme to 
solve the Einstein's equations, energy-momentum conservation equations and baryonic 
conservation equation. The code is similar to the one originally developed by May and 
White (1967). Here we discuss the capabilities of the code that make it well suited for 
numerical relativity on a personal computer and some caveats based on the experiments 
we have made with it. 



1. Introduction 

In this work we have reproduced the code of May and White (1967) (hereafter 
MW) using the same numerical techinques (with some little modifications) known 
from their pioneering work. We had also tested some modifications suggested by 
Van Riper (1979) (hereafter VR). The output of our code can be used as an initial 
Cauchy data for more sophisticated relativistic codes (see for example Baumgarte 
et al. 1995). The code developed here altough based on old numerical techniques is 
simple to implement and little time consuming, which make it best suited for general 
relativistic simulations on a single personal computer. Here we briefly resume the 
numerical caveats on a simulation of the hydrodynamic collapse of a 21 solar masses 
star that forms a black hole. The complete results of astrophysical interest will be 
published elsewhere (Ghezzi & Letelier 2003). 

2. Numerical methods and initial conditions 

We initially assumed uniform density perfect fluid balls, with adiabatic equation 
of state P = Kp" 1 , where 7 is the adiabatic coefficient which is set equal to 5/3. 
It is also possivel to implement other equations of state in our code. The initial 
configuration resembles a newtonian star or a newtonian star-core collapsing under 
the gravity when nuclear fuel is exhausted. During the time evolution a strong field 
regime is achieved during which a neutron star or a black hole (a trapped surface) 



7, 2008 2:42 WSPC/Trim Size: 9.75in x 6.5in for Proceedings 



main 



2 

may be formed depending on the initial conditions. The formation of the black 
hole was checked by the collapse of the gtt metric coefficient (collapse of the "lapse 
function"). Alternatively, the black hole formation could be checked by observing 
the convergence of the ligth ray paths that pass near a Lagrangian coordinate /z, 
such that 2GM(fi)/Rc 2 > 1. 

It is possivel to assume arbitrary density or energy profiles as initial data sets, 
altough we must observe that could be troublesome if discontinuities are present on 
the initial data. 

We have used a Lax-Wendroff scheme (two step Ritchmeyer's version) coupled 
with a Crank-Nicolson (1947) algorithm to solve the system of equations (similar 
to the techniques described in MW, and in VR). We have made tests using the 
viscocities given in VR and in MW. The boundary conditions are that the pressure 
is zero at surface, the mass of the star is zero at the centre of the star, and the 
coefficient g u of the metric is equal to one at the surface of the star. This last 
boundary condition assures that the coordinate time are synchronised with the 
clock of a comoving observer at the surface of the collapsing star. The results 
obtained from one run test are shown in the figures (a)-(d). 

A novelty feature of our code is that it don't break down when a black hole 
forms. This is because the code was built to self-adjust the time step in order to 
get the desired accuracy, and eventually the time step becames exceedingly small 
when a black hole is obtained as the final outcome of the simulation. The code is 
very fast and accurate. 

3. Discussion and Conclusion 

Using the VR viscocity or the MW viscocity seems to give no different results using 
a low number of radial zones (~ 70). However, when a larger number of zones 
is set the MW viscocity seems to work better if it is multiplied by some constant 
factor greater than one. This is necessary to enlarger the shock waves through the 
zones, because when a large number of radial zones is used the shocks have a steeper 
representation, and this in turns could lead to numeric dispersion and consequently 
to a low accuracy. 

It is possible to make very fast simulations, in the present case for example, 
using 200 radial zones the code take <~ 10 — 30 minutes to run on a Pentium IV, 1.7 
Mhz. Altough, to go further in the time evolution of the collapse could take much 
more time. The total accuracy at the end of the simulation is ~ 99.5 %. However, 
depending mostly on the choice of the viscocity parameter the accuracy could be 
worst. 

We must point out that the code only depart from excellent energy conservation 
only when a trapped surface begins to form, so we think the precision obtained is 
quit satisfactory. 
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(a) (b) 

Figure 1. The figures displays a time sequence of snapshots as a function of the Lagrangian 
coordinate [i for the density (Fig. a), and for the metric coefficient gtt (Fig. b). 




(i [solar masses] \L [solar massesj 



(c) (d) 

Figure 2. Snapshots for the total mass or energy (Fig. c), and for the radius of the configuration 
(Fig. d). 
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